#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Implement several metrics for subspaces."""

__revision__ = "$Revision$"

import numpy as np
import math

from .projections import calcProjectionMatrix


__all__ = ["calc_principal_angles",
           "calc_chordal_distance_from_principal_angles",
           "calc_chordal_distance", "calc_chordal_distance_2"]


# TODO: I think calc_principal_angles is not correct when matrix1 e matrix2
# have different sizes. At least obtaining the chordal distance from the
# principal angles does not work when matrix1 and matrix2 have different
# shapes.
def calc_principal_angles(matrix1, matrix2):
    """Calculates the principal angles between `matrix1` and `matrix2`.

    Parameters
    ----------
    matrix1 : 2D numpy array
        A 2D numpy array.
    matrix2 : 2D numpy array
        A 2D numpy array.

    Returns
    -------
    princ_angles : 1D numpy array
        The principal angles between `matrix1` and `matrix2`.

    See also
    --------
    calc_chordal_distance_from_principal_angles

    Examples
    --------
    >>> A = np.array([[1, 2], [3, 4], [5, 6]])
    >>> B = np.array([[1, 5], [3, 7], [5, -1]])
    >>> print(calc_principal_angles(A, B))
    [ 0.          0.54312217]
    """
    # First we need to find the orthogonal basis for matrix1 and
    # matrix2. This can be done with the QR decomposition. Note that if
    # matrix1 has 'n' columns then its orthogonal basis is given by the
    # first 'n' columns of the 'Q' matrix from its QR decomposition.
    Q1 = np.linalg.qr(matrix1)[0]
    Q2 = np.linalg.qr(matrix2)[0]

    # TODO: Test who has more columns. Q1 must have dimension grater than
    # or equal to Q2 so that the SVD can be calculated in the order below.
    #
    # See the algorithm in
    # http://sensblogs.wordpress.com/2011/09/07/matlab-codes-for-principal-angles-also-termed-as-canonical-correlation-between-any-arbitrary-subspaces-redirected-from-jen-mei-changs-dissertation/
    S = np.linalg.svd(
        Q1.conjugate().transpose().dot(Q2), full_matrices=False)[1]

    # The singular values of S vary between 0 and 1, but due to
    # computational impressions there can be some value above 1 (by a very
    # small value). Below we change values greater then 1 to be equal to 1
    # to avoid problems with the arc-cos call later.
    S[S > 1] = 1  # Change values greater then 1 to 1

    # The singular values in the matrix S are equal to the cosine of the
    # principal angles. We can calcuate the arc-cosine of each element
    # then.
    return np.arccos(S)


def calc_chordal_distance_from_principal_angles(principalAngles):
    """Calculates the chordal distance from the principal angles.

    It is given by the square root of the sum of the squares of the sin of
    the principal angles.

    Parameters
    ----------
    principalAngles : 1D numpy array
        Numpy array with the principal angles.

    Returns
    -------
    chord_dist : int
        The chordal distance.

    See also
    --------
    calc_principal_angles,
    calc_chordal_distance,
    calc_chordal_distance_2

    Examples
    --------
    >>> A = np.arange(1, 9.)
    >>> A.shape = (4, 2)
    >>> B = np.array([[1.2, 2.1], [2.9, 4.3], [5.2, 6.1], [6.8, 8.1]])
    >>> princ_angles = calc_principal_angles(A, B)
    >>> print(calc_chordal_distance_from_principal_angles(princ_angles))
    0.473867859572
    """
    return np.sqrt(np.sum(np.sin(principalAngles) ** 2))


def calc_chordal_distance(matrix1, matrix2):
    """Calculates the chordal distance between the two matrices

    Parameters
    ----------
    matrix1 : 2D numpy array
        A 2D numpy array.
    matrix2 : 2D numpy array
        A 2D numpy array.

    Returns
    -------
    chord_dist : int
        The chordal distance.

    Notes
    -----
    Same as :func:`calc_chordal_distance_2`, but implemented differently.

    See also
    --------
    calc_chordal_distance_2,
    calc_chordal_distance_from_principal_angles

    Examples
    --------
    >>> A = np.arange(1, 9.)
    >>> A.shape = (4, 2)
    >>> B = np.array([[1.2, 2.1], [2.9, 4.3], [5.2, 6.1], [6.8, 8.1]])
    >>> print(calc_chordal_distance(A, B))
    0.473867859572
    """
    Q1 = np.linalg.qr(matrix1)[0]
    Q2 = np.linalg.qr(matrix2)[0]

    #ncols = matrix1.shape[1]  # Deve ser igual a matrix2.shape[1].

    # As primeiras ncols colunas de Q1 e Q2 formam a base ortonormal de
    # ran(matrix1) e ran(matrix2), respectivamente
    Q1_sqr = Q1.dot(Q1.conjugate().transpose())
    Q2_sqr = Q2.dot(Q2.conjugate().transpose())
    return np.linalg.norm(Q1_sqr - Q2_sqr, 'fro') / math.sqrt(2.)


def calc_chordal_distance_2(matrix1, matrix2):
    """Calculates the chordal distance between the two matrices

    Parameters
    ----------
    matrix1 : 2D numpy array
        A 2D numpy array.
    matrix2 : 2D numpy array
        A 2D numpy array.

    Returns
    -------
    chord_dist : int
        The chordal distance.

    Notes
    -----
    Same as :func:`calc_chordal_distance`, but implemented differently.

    See also
    --------
    calc_chordal_distance,
    calc_chordal_distance_from_principal_angles

    Examples
    --------
    >>> A = np.arange(1, 9.)
    >>> A.shape = (4, 2)
    >>> B = np.array([[1.2, 2.1], [2.9, 4.3], [5.2, 6.1], [6.8, 8.1]])
    >>> print(calc_chordal_distance_2(A, B))
    0.473867859572
    """
    return (np.linalg.norm(calcProjectionMatrix(matrix1)
                           - calcProjectionMatrix(matrix2), 'fro')
            / math.sqrt(2))
